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Abstract — We describe the construction of a spherical wavelet 
analysis through the inverse stereographic projection of the 
Euclidean planar wavelet framework, introduced originally by 
Antoine and Vandergheynst and developed further by Wiaux 
et al. Fast algorithms for performing the directional continuous 
wavelet analysis on the unit sphere are presented. The fast 
directional algorithm, based on the fast spherical convolution 
algorith m de veloped by Wandelt and Gorski, provides a saving 
of G(y/N P i X ) over a direct quadrature implementation for Ar pix 
pixels on the sphere, and allows one to perform a directional 
spherical wavelet analysis of a 10 6 pixel map on a personal 
computer. 

Index Terms — Wavelet transforms, spheres, convolution. 



I. Introduction 

WAVELET analysis has proven useful in many ap- 
plications due to the ability of wavelets to resolve 
localised signal content in both scale and space. Many of 
these applications, however, are restricted to data defined in 
Euclidean space: the 1 -dimensional line, the 2-dimensional 
plane and, occasionally, higher dimensions. Nevertheless, data 
are often measured or defined on other manifolds, such as the 
2- sphere. Examples where data are measured on the sphere 
are found in astrophysics (e.g. [1], [2]), planetary science 
(e.g. [3]-[5]), geophysics (e.g. [6]-[8]), computer vision (e.g. 
[9]) and quantum chemistry (e.g. [10], [11]). To realise the 
potential benefits that may be provided by wavelets in such 
settings, ordinary Euclidean wavelet analysis must be extended 
to spherical geometry. 

A number of attempts have been made to extend wavelets 
to the unit sphere. Discrete second generation wavelets on the 
sphere that are based on a multiresolution analysis have been 
developed [12], [13]. Other authors have focused on the con- 
tinuous wavelet transform on the sphere. A number of works 
construct a solution using a harmonic approach [14]— [17], 
however these solutions suffer from the poor localisation of the 
spherical harmonic functions. Others adopt a tangent bundle 
viewpoint [18], [19], thereby avoiding the necessity to define 
a dilation operator on the sphere. A satisfactory extension of 
the continuous wavelet transform to the sphere is defined by 
[20], however this construction requires an abstract dilation 
parameter that must satisfy a number of ad hoc assumptions. 
More recently, a consistent and satisfactory framework for 
wavelets defined on the unit sphere has been constructed and 
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developed by [21]-[27]. Moreover, this construction is derived 
entirely from group theoretic principles and inherently satisfies 
a number of natural requirements. We consider the continuous 
spherical wavelet transform (CSWT) developed in these last 
works. For a more detailed review of the attempts made at 
constructing a wavelet transform on the unit sphere see [21], 
[22], [28]. 

Current and future data- sets defined on the sphere are 
of considerable size. The current Wilkinson Microwave 
Anisotropy Probe (WMAP) data of the cosmic microwave 
background (CMB) contain approximately 3 x 10 6 pixels on 
the sphere, whereas the forthcoming Planck CMB mission 
will generate maps with approximately 50 x 10 6 pixels. 
Fast algorithms are therefore required to perform the CSWT 
on practical data- sets. A semi-fast algorithm to implement 
the CSWT is presented in [23] and implemented in the 
YAWTb 1 (Yet- Another- Wavelet-Toolbox) Matlab wavelet tool- 
box (which also makes use of the SpharmonicKit 2 ). However, 
this implementation is restricted to an equi-angular tessellation 
of the sphere. The beauty of this tessellation is its simplicity 
and ability to be easily represented in matrix form. However, 
the pixels of an equi-angular tessellation are densely spaced 
about the poles and do not have equal areas. Other tessellations 
of the sphere also exist, such as those constructed to minimise 
some energy measure [29]-[31] or those constructed for more 
practical or numerical purposes (for example the IGLOO 3 [32], 
HEALPix 4 [33] and GLESP 5 [34] schemes). There is thus a 
need for a fast implementation of the CSWT that is not tied 
to any particular tessellation of the sphere. 

We fill this void by presenting a fast algorithm for perform- 
ing the directional CSWT. The CSWT at a particular scale is 
essentially a spherical convolution, hence we may apply the 
fast spherical convolution algorithm developed by [35] to eval- 
uate the wavelet transform. The algorithm is posed in harmonic 
space and thus is independent of the underlying tessellation of 
the sphere, (although an iso-latitude tessellation does enable 
faster spherical harmonic transforms, thereby increasing the 
speed of the algorithm). The framework supports both non- 
azimuthally symmetric spherical wavelets 6 and a decomposi- 
tion that employs anisotropic dilations, however no synthesis 
is possible when anisotropic dilations are incorporated. For 



^ttp : / /www . f yma . ucl . ac . be/pro ject s/yawtb/ 
2 http: //www. cs .dartmouth.edu/~geelong/sphere/ 
3 http : //www .mrao . cam. ac . uk/pro ject s/cpac/ igloo/ 
4 http : / /healpix . jpl . nasa . gov/ 
5 http : / /www .glesp.nbi.dk/ 

6 Azimuthally symmetric spherical wavelets are also often referred to as 



axisymmetric wavelets. 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. -, NO. -, JUNE 2005 



an illustration of a spherical wavelet analysis on a practical 
problem of considerable size we refer the reader to our recent 
works to test the CMB for deviations from Gaussianity [36], 
[37] and to detect the integrated Sachs-Wolfe (ISW) effect 
[38]. Both of these works involve performing 1000 Monte 
Carlo simulations and would not have been feasible without a 
fast directional CSWT algorithm. 

The remainder of this paper is structured as follows. In 
section II we describe the CSWT in the framework presented 
by [27]. For a more complete treatment of the spherical 
wavelet transform in this framework and the correspondence 
between spherical and Euclidean wavelets we recommend that 
the reader refer to [27]. We also present an extension to 
anisotropic dilations, however in this case the basis functions 
are not strictly wavelets hence perfect reconstruction is not 
possible. Various algorithms to perform the CSWT are de- 
scribed and then compared in section III. An application of 
our implementation is demonstrated in section IV. Concluding 
remarks are made in section V. 

II. The continuous spherical wavelet transform 

To extend Euclidean wavelet analysis to spherical geometry 
a number of requirements must be satisfied: (i) the signals 
and wavelets must live fully on the unit sphere; (ii) the 
transform must involve local dilations of some kind on the 
unit sphere; and (iii) the spherical wavelet transform should 
reduce locally to the Euclidean transform on the tangent plane 
(i.e. the Euclidean limit must be satisfied) [20]. The final 
requirement is intuitively obvious; the sphere is asymptotically 
flat, hence any spherical wavelet transform should match the 
planar Euclidean transform on small scales, or equivalently, 
for a large radius of curvature. 

The spherical wavelet transform developed in [21]-[26] 
satisfies all of these requirements, moreover each requirement 
naturally follows from the construction. The construction 
of this transform is derived entirely from group theoretical 
principles. However, in a recent work by [27] this formalism 
is reintroduced independently of the original group theoretic 
formalism, in an equivalent, practical and self-consistent ap- 
proach. We adopt this approach herein. 

The correspondence principle between spherical and Eu- 
clidean wavelets is developed by [27], relating the concepts 
of planar Euclidean wavelets to spherical wavelets through a 
stereographic projection. We use the stereographic projection 
to define affine transformations on the unit sphere that facil- 
itate the construction of a wavelet basis on the unit sphere. 
The spherical wavelet transform may then be defined as the 
projection on to this basis, where the spherical wavelets must 
satisfy the appropriate admissibility criterion to ensure perfect 
reconstruction. 

A. Stereographic projection 

In order to construct a correspondence between wavelets on 
the plane (R 2 ) and sphere (S 2 ) a projection operator between 
the two spaces must be chosen. It is shown in [27] that 
the stereographic projection is the unique unitary, radial and 
conformal diffeomorphism between the sphere and the plane. 




South pole 



Fig. 1. Stereographic projection of the sphere onto the plane. 

The stereographic projection is defined by projecting a point 
on the unit sphere to a point on the tangent plane at the north 
pole, by casting a ray though the point and the south pole. 
The point on the unit sphere is mapped on to the intersection 
of this ray and the tangent plane (see Fig. 1). Formally, we 
may define the stereographic projection operator as II : uj — > 
x = = (r(<9), (j)) where r = 2 tan(0/2), u = (0, <j>) G S 2 
denotes spherical coordinates with colatitude and longitude 
(j) and x G R 2 is a point in the plane, denoted here by the 
polar coordinates (r, 4>). The inverse operator is II -1 : x — > 
u = n _1 x = (0(r),0), where 0(r) = 2 tan- 1 (r/2). 

Following the formulation of [27] again, we define the 
action of the stereographic projection operator on functions on 
the plane and sphere. We consider the space of square inte- 
grable functions in L 2 (R 2 , d 2 x) on the plane and L 2 (S 2 , dfl) 
on the unit sphere, where dft = sinOdOdcj) is the usual 
rotation invariant measure on the sphere. The action of the 
stereographic projection operator II : s G L 2 (S 2 , dQ) — > p = 
lis G L 2 (R 2 , d 2 x) on functions is defined as 

p(r, 0) = (IL>)(r, 0) = (1 + r 2 /4)- 1 s(0(r), 0) . (1) 

The inverse stereographic projection operator II -1 : p G 
L 2 (R 2 , d 2 x) -> s = Ii~ l p G L 2 (S 2 , dfi) on functions is 
then 

s(0, <j>) = (n- 1 p)(^, </») = [! + tan 2 (0/2)Mr(0), <f>) . (2) 

The pre-factors introduced ensure that the L 2 -norm of func- 
tions through the forward and inverse projections are con- 
served. In the Euclidean limit, the stereographic projection 
and inverse naturally reduce to the identity operator [21]. 

B. Affine transformations on the sphere 

A wavelet basis is constructed on the unit sphere in sec- 
tion II-C by applying the spherical extension of Euclidean 
translations and dilations to mother wavelets defined on the 
unit sphere. The extension of these affine transformations to 
the sphere, facilitated by the stereographic projection operator, 
are defined here. 

The natural extension of Euclidean translations on the unit 
sphere are rotations. These are characterised by the elements 
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of the rotation group SO (3), which we parameterise in terms 
of the three Euler angles p = (a,/3, 7). 7 The rotation of a 
square-integrable function s G L 2 (S 2 , dfi) is defined by 



(3) 



Dilations on the unit sphere are constructed by first lift- 
ing the sphere to the plane by the stereographic projection, 
followed by the usual Euclidean dilation in the plane, be- 
fore re-projecting back onto the sphere. We generalise to 
anisotropic dilations on the sphere (a similar anisotropic 
dilation operator on the sphere has been independently pro- 
posed by [39]), however in this setting we do not achieve 
a wavelet basis and hence cannot synthesise our original 
signal. We define the anisotropic Euclidean dilation operator 
in L 2 (R 2 , d 2 x) as 



[d(a, b)p] y) = ar^b- 1 ' 2 p^x, b^y) , 



(4) 



for the non-zero positive scales a, b G R+. The a~ 1 / 2 b~ 1 / 2 
normalisation factor ensures the L 2 -norm is preserved. The 
spherical dilation operator £>(a, b) : s(0, <p) — > [V(a, b)s] (0, (j)) 
in L 2 (£ 2 , dfl) is defined as the conjugation by II of the 
Euclidean dilation d(a, b) in L 2 (R 2 , d 2 x) on the tangent plane 
at the north pole: 



V(a,b) = n _1 d(a,6)n. 



(5) 



The norm of functions in L 2 (S 2 , dQ) is preserved by the 
spherical dilation as both the stereographic projection oper- 
ator and Euclidean dilations preserve the norm of functions. 
Extending the isotropic spherical dilation operator defined by 
[27] to anisotropic dilations, we obtain 

[V{a,b)s](u) = [A(a,6,0,0)] 1/2 *(ui/a,i/b) , (6) 
where cj a ,b = (0 a ,b, <t>a,b), 



tan((9 a , 6 /2) = tan((9/2) 



b 2 sin 2 



and 



tan(0 o b ) = - tan(0) 
a 



For the case where a = 6 the anisotropic dilation reduces to the 
usual isotropic case defined by tan(0 a /2) = atan(0/2) and 
<j) a = cj). The A (a, 6, 0, 0) cocycle term follows from the factors 
introduced in the stereographic projection of functions to pre- 
serve the L 2 -norm. Alternatively, the cocycle may be derived 
explicitly to preserve the L 2 -norm when the stereographic 
projection of functions do not have these pre-f actor terms. The 
cocycle of an anisotropic spherical dilation is defined by 



where 



4a 3 b 3 

X(a. 6, 0, 6) = o 

V J (A_cos0 + A+) 2 



A+ = a 2 b 2 ± a 2 sin 2 6 ± b 2 cos 2 . 



(7) 



7 We adopt the zyz Euler convention corresponding to the rotation of a 
physical body in a fixed co-ordinate system about the z, y and z axes by 7, 
(3 and ot respectively. 



For the case where a = b the anisotropic cocycle reduces to 
the usual isotropic cocycle 

4a 2 

A(a ' a ' = [(a 2 -l)cos^ + a 2 + l] 2 * 

Although the ability to perform anisotropic dilations is of 
practical use, we do not achieve a wavelet basis in this setting. 
In the anisotropic setting a bounded admissibility integral 
cannot be determined (even in the plane), thus the synthesis of 
a signal from its coefficients cannot be performed. This results 
from there being no direct means of evaluating the proper 
measure in the absence of a group structure. The projection of 
a signal onto basis functions undergoing anisotropic dilations 
may be performed in an analogous manner to the following 
discussion of the wavelet transform. However, since these basis 
functions are not wavelets we restrict the following discussion 
to isotropic dilations. 

C. Wavelet transform 

A wavelet basis on the unit sphere may now be constructed 
from rotations and isotropic dilations (where a = b) of a 
mother spherical wavelet ijj G L 2 (S 2 , dfl). The corresponding 
wavelet family {^ a , P = R(p)V(a,a)ip, p G SO(3), a G R+} 
provides an over-complete set of functions in L 2 (S 2 , dQ). The 
CSWT of s G L 2 (S 2 , dfl) is given by the projection onto each 
wavelet basis function in the usual manner, 



W^p) 



(8) 



where the * denotes complex conjugation. 

The transform is general in the sense that all orientations 
in the rotation group SO (3) are considered, thus directional 
structure is naturally incorporated. It is important to note, 
however, that only local directions make any sense on £ 2 . 
There is no global way of defining directions on the sphere 8 - 
there will always be some singular point where the definition 
fails. 

The synthesis of a signal on the unit sphere from its wavelet 
coefficients is given by 

s(u) = I dp r ^ W$(a,p) [R(p)L^ a ](u) , (9) 

JSO(3) JO a 

where dp = sin/3dad^d7. The operator in L 2 (S 2 , dfl) 
is defined by the action 



(10) 



on the spherical harmonic coefficients of functions 
/ G L 2 (S 2 , dfl), where is defined below. The hat 
denotes the spherical harmonic coefficients 

hm= I doy/ m H/H = (y, m |/) 

Js 2 

of the decomposition 

00 i 

£=0 m =-£ 

8 There is no differentiable vector field of constant norm on the sphere and 
hence no global way of defining directions. 
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We adopt the Condon- Shortley phase convention where the 
normalised spherical harmonics are defined by 



where P™(x) are the associated Legendre functions. Using 
this normalisation the orthogonality of the spherical harmonic 
functions is given by 



/ 

Js 2 



dSlY£ m (u;)Y^ mf (u;) 8w&mm' , 



(11) 



where Sij is Kronecker delta function. In order to ensure the 
perfect reconstruction of a signal synthesised from its wavelet 
coefficients, one requires the admissibility condition 



8tt 5 



T E / -T I ?< 00 (12) 



to hold for all £ G N. A proof of the admissibility criterion 
is given by [27]. Practically, it is difficult to apply (12) 
directly, thus a necessary (and almost sufficient) condition for 
admissibility is the zero-mean condition [21] 



S 2 1 + cos 







(13) 



In the Euclidean limit this condition naturally reduces to the 
necessary zero-mean condition for Euclidean wavelets [27]. 

D. Correspondence principle and spherical wavelets 

The correspondence principle between spherical and Eu- 
clidean wavelets states that the inverse stereographic projection 
of an admissible wavelet on the plane yields an admissible 
wavelet on the unit sphere. This result is proved by [27]. 
Hence, mother spherical wavelets may be constructed from 
the projection of mother Euclidean wavelets on the plane: 

^M = [n-V^]M, (14) 

where ij^2 G L 2 (R 2 , d 2 x) is an admissible wavelet in the 
plane. Directional spherical wavelets may be naturally con- 
structed in this setting - they are simply the projection of 
directional Euclidean planar wavelets on to the sphere. 

We give examples of three spherical wavelets: the spherical 
Mexican hat wavelet (SMHW); the spherical butterfly wavelet 
(SBW); and the spherical real Morlet wavelet (SMW). These 
spherical wavelets are illustrated in Fig. 2. Each spherical 
wavelet is constructed by the stereographic projection of the 
corresponding Euclidean wavelet onto the sphere, where the 
Euclidean planar wavelets are defined by 



SMHW 



(r, 



-(2-r 2 )e" r2 / 2 



iif W (x,y)=xe 



-(x 2 +y 2 )/2 



and 



(x;k) =Re(e lk - x /^e-W 2 / 2 ) 



respectively, where k is the wave vector of the SMW. The 
SMHW is proportional to the Laplacian of a Gaussian, 
whereas the SBW is proportional to the first partial derivative 



of a Gaussian in the ^-direction. The SMW is a Gaussian 
modulated sinusoid, or Gabor wavelet. 

A full directional wavelet analysis on the unit sphere for 
large data sets has previously been prohibited by the computa- 
tional infeasibility of any implementation. The computational 
burden of computing many orientations may be reduced by 
using steerable wavelets, for which any continuous orientation 
can be computed from a small number of basis orientations 
[27]. This is achieved since steerable wavelets have a limited 
azimuthal band limit and may thus be represented as a finite 
sum of trigonometric exponentials [27]. However, in this case 
one must still compute the initial transform for more than one 
orientation, so although the computational burden is reduced, 
it is still significant. Moreover, we also require a fast approach 
for general non-steerable, directional wavelets. We address this 
problem in the following section by presenting a fast algorithm 
to perform the directional CSWT. 

III. Algorithms 

A range of algorithms of varying computational efficiency 
and numerical accuracy are presented to perform the CSWT 
described in section II. We implement these algorithms in For- 
tran 90 and subsequently compare computational complexity 
and typical execution time. The synthesis of a signal from its 
wavelet coefficients is not considered any further. Without loss 
of generality we consider only a single dilation (i.e. fixed a 
and b). 

A. Tessellation schemes 

It is necessary to discretise both the spherical coordinates 
of a function defined on the unit sphere and also the Euler 
angle representation of the SO (3) rotation group. The fast 
algorithms we present are performed in harmonic space and 
hence are tessellation independent, provided an appropriate 
spherical harmonic transform is defined for the tessellation. 
However, the semi-fast algorithm is restricted to an equi- 
angular tessellation of the sphere. The various tessellation 
schemes adopted are defined below. 

The equi-angular tessellation (also known as the equidistant 
cylindrical projection (ECP)) of the spherical coordinates is 



: < n e < N e 



defined by C = {0 ne = ^f^ ne - 
1, < < — 1}. Let A^pi x = NqN^ denote the number 
of pixels in the tessellation. 

We also consider the HEALPix tessellation scheme since it 
is commonly used for astrophysical data-sets of the CMB. 
Pixels in the HEALPix scheme are of equal area and are 
located on rings of constant latitude (the latter feature enables 
fast spherical harmonic transforms on the pixelised grid). We 
refer the reader to [33] for details of the HEALPix scheme and 
here just define the HEALPix grid in terms of pixel indices: 
TL = {(<9, <j>) p = (6 P , </> p ) : < p < Np ix - 1}. The HEALPix 
resolution is parameterised by iV s ide, where 7V p i x = 127V S id e 2 . 
It should be noted that an exact quadrature formula does not 
exist for the HEALPix tessellation, thus spherical harmonic 
transforms are necessarily approximate. This is not the case 
for the ECP or other practical tessellations (e.g. IGLOO and 
GLESP) where exact quadrature may be performed. 
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(a) Spherical Mexican hat wavelet (b) Spherical butterfly wavelet (SBW) (c) Spherical real Morlet wavelet (SMW) 

(SMHW) 

Fig. 2. Spherical wavelets at scale a = b = 0.2. Wavelet maps are displayed in the Mollweide projection, where the wavelets have been rotated down from 
the north pole for ease of observation. The SMW is plotted for wave vector k = (10, 0) T . 



The Euler angle domain of the spherical wavelet coefficients 
is in general arbitrary, however we use the equi-angular dis- 
cretisation defined by £ 1 {a na p np = ^f, ^ = 



-jT* : < n a < N a - 1, < np < Np - 1, < n 7 < N 7 - 
1}. Our fast algorithm, however, requires (for convenience) the 



tessellation £2 = {ct na = 



27rn Q 



5 fin 



2itnQ 
~N~n 



7n 7 



2irn 1 



< n a < N a - 1, < np < 2Np - 1, < n 7 < N 1 - 1}, 
where the fi sampling is repeated. Evaluating (3 over the range 
to 2tt is redundant, covering the SO (3) manifold exactly 
twice. Nonetheless, the use of our fast algorithm requires this 
range. Such an approach is not uncommon, as [40], [41] also 
oversample a function of the sphere in the direction in 
order to develop fast spherical harmonic transforms on equi- 
angular grids. The overhead associated with our inefficient 
discretisation is more than offset by the fast algorithm it 
affords, as described in section III-E. 

B. Direct case 

The CSWT defined by (8) may be implemented directly 
by applying an appropriate quadrature rule. Using index sub- 
scripts to denote sampled signals, the direct CSWT implemen- 
tation is given by 



ip )riu,np,n 1 
iVpix-1 



E 

p=0 



R 



2im a nnp 27rn 7 



s p w p , (15) 



where the pixel sum is over all the pixels of any chosen 
tessellation. The weights for the ECP C grid are given by 

whereas the equal pixel areas of 



2tt 2 sin 6 n 



the HEALPix H grid ensure the pixel weights, given by 



— , are independent of position. 



Discretisation techniques other than the plain Riemann sum 
used in (15) would be beneficial only if additional regularity 
conditions are imposed on the signal s [23]. It is also possible 
to choose other weights to achieve a better approximation of 
the integral. An example of a different equi-angular discreti- 
sation and a different choice for the weights is given by the 
sampling theorem for band-limited functions on the sphere 
developed by [41]. 

Evaluation of (15) requires the computation of a 
2-dimensional summation, evaluated over a 3 -dimensional 



grid. We assume the number of samples for each discretised 
angle, except 7, is of the same order N. Typically the number 
of samples in the 7 direction is of a much lower order, so 
we treat this term separately. The complexity of the direct 
algorithm is 0(N 1 N A ). 



C. Semi-fast case 

We rederive here the semi-fast implementation of the CSWT 
described by [23] and implemented using the YAWTb Matlab 
wavelet toolbox and SpharmonicKit. This algorithm involves 
performing a separation of variables so that one rotation may 
be performed in Fourier space. The algorithm is restricted to 
the equi-angular grid C (in essence pixels must only be defined 
on equal latitude rings, however some form of interpolation 
and down- sampling is then required to extract samples for 
equal longitudes). 

Firstly, the a rotation is represented by shifting the corre- 
sponding wavelet samples to give 



(W^)n a ,n f3 ,n 1 

N e -l N4-I r 

E E 

n g =0 n,f,=0 



R\0 



Tmp_ 27rn 7 



$n e ,n6 W ne 1 



where the index is extended periodically with period N^. 
The discrete space convolution theorem may then be applied to 
represent the inner summation as the inverse discrete Fourier 
transform (DFT) of the product of the wavelet and signal DFT 
samples (note that only a 1 -dimensional DFT is performed in 
the azimuthal direction): 



ijj )n a ,np,n 1 



Ne-l 



n =O ^ 

X F(s) 



1 

n e ,k 



N^ 



E^* 



k=0 



R[0, 



imp 27m n 



n e ,k 



(16) 



where ^"(-) nj fe denotes 1 -dimensional DFT coefficients. A 
fast Fourier transform (FFT) may then be applied to evaluate 
simultaneously all of the n a terms of the expression enclosed 
in the curly braces in (16). A final summation (integral) 
over ng produces the spherical wavelet coefficients for a 
given np and n 7 , for all n a . Applying an FFT to evaluate 
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simultaneously one summation rapidly, reduces the complexity 
of the CSWT implementation to 0(N 7 N 3 log 2 N). 

D. Fast azimuthally symmetric case 

The fast azimuthally symmetric CSWT algorithm is posed 
in harmonic space, where sg m = (Y^ m | s) are the spherical 
harmonic coefficients of a function s G L 2 (S 2 , d£7), as 
described in section II-C. For the special case where the 
wavelet is azimuthally symmetric (i.e. invariant under az- 
imuthal rotations), it is essentially only a function of and 
may be represented in terms of its Legendre expansion. In this 
case the harmonic representation of the wavelet coefficients 
is given by the product of the signal and wavelet spherical 
harmonic coefficients: 



4tt 



2£- 



^£0 Sin 



(17) 



noting that the harmonic coefficients of an azimuthally sym- 
metric wavelet are zero for m ^ 0. In practice, one requires 
that at least one of the signals, usually the wavelet, has a 
finite band limit so that negligible power is present in those 
coefficients above a certain ^ max . We then only need to 
consider £ < £ max (a detailed discussion of the determination 
of tax is presented in [25]). Once the spherical harmonic 
representation of the wavelet coefficients is calculated by (17), 
the inverse spherical harmonic transform is applied to compute 
the wavelet coefficients in the Euler domain. The complexity 
of the fast isotropic CSWT algorithm is dominated by the 
spherical harmonic transforms. For a tessellation containing 
pixels on rings of constant latitude, a fast spherical harmonic 
transform may be performed (see e.g. [33], [42], [43]). This 
reduces the complexity of the spherical harmonic transform 
from G(N 4 ) to 0(N 3 ) = 0(N^ 3/2 ). For certain tessellation 
schemes fast spherical harmonic transforms of lower complex- 
ity are also available, however these are related directly to 
the tessellation (e.g. [40]-[43]). In particular, the algorithm 
developed by [40], [41] for the ECP tessellation scales as 
G(N 2 log TV). 

The fast azimuthally symmetric CSWT algorithm is posed 
purely in harmonic space and consequently the algorithm 
is tessellation independent. However, we are restricted to 
azimuthally symmetric wavelets and lose the ability to perform 
the directional analysis inherent in the wavelet transform 
construction. 

E. Fast directional case 

We present the most general fast directional CSWT algo- 
rithm for non-azimuthally symmetric wavelets, i.e. steerable 
and directional wavelets, in this section. Again, the algorithm 
is posed purely in harmonic space and so is tessellation 
independent. We do, however, use the equi-angular £ 2 dis- 
cretisation of the wavelet coefficient domain, although other 
discretisations may be used if FFTs are also defined on these 
grids. The CSWT at a particular scale (i.e. a particular a 
and b) is essentially a spherical convolution, hence we apply 
the fast spherical convolution algorithm proposed by [35] to 
evaluate the wavelet transform. The algorithm proceeds by 



factoring the rotation into two separate rotations, each of which 
involves only a constant polar rotation component. Azimuthal 
rotations may then be performed in harmonic space at far less 
computation expense than polar rotations. We subsequently 
rederive the fast spherical convolution algorithm developed by 
[35], as applied to our application of evaluating the CSWT. 
The harmonic representation of the CSWT is first presented, 
followed by the discretisation and fast implementation. 

1 ) Harmonic formulation: Substituting the spherical har- 
monic expansions of the wavelet and signal into the wavelet 
transform defined by (8), and noting the orthogonality of the 
spherical harmonics described by (11), yields the harmonic 
representation 



E E E 



£=0 m=-£ 



DLm'faPil) fan 



S£m • (18) 



Again, we assume negligible power above £ max in at least 
one of the signals, usually the wavelet, so that the outer 
summation is truncated to -£ max . The additional summation 
over m' and the £>^ m / Wigner rotation matrices that are 
introduced characterise the rotation of a spherical harmonic, 
noting that a rotated spherical harmonic may be represented 
by a sum of harmonics of the same £ [11], [44]: 



[R(a,(3 n )Y £m ](u;)= Di, m (a,f3^)Y£ m ,(u). 

m' = -£ 

The Wigner rotation matrices may be decomposed as [11], 
[44] 



(19) 



where the real polar d-matrix is defined by [44] 

min(£-\-m,£— m') 

<L m >(P)= E ("!)' 

£ = m ax ( , rn — rn ' ) 

[(£ + m)\ (£ - m)\ (£ + m')\ (£ - m')\ ] 1/2 
X (i + m - t)\ (£ - m' - t)\ (t + m' - m)\ t\ 



cos 



-i 2£+m-m' ' -It r 



sin 



1 m f —m-\-2t 



and the sum over t is defined so that the arguments of the 
factorials are non-negative. Recursion formulae are available 
to compute rapidly the Wigner rotation matrices in the basis 
of either complex [10], [45] or real [46], [47] spherical 
harmonics. We employ the recursion formulae described in 
[45] in our implementation. The decomposition shown in (19) 
is exploited by factoring the rotation R(a,/3,j) into two 
separate rotations, both of which only contain a constant ±7r/2 
polar rotation: 

R(a,f3,j) = R(a-n/2, -tt/2, 0) 

xR(0, tt/2, 7 + 7T/2). (20) 



MCEWEN et al: FAST DIRECTIONAL CSWT ALGORITHMS 



7 



By factoring the rotation in this manner and applying the 
decomposition described by (19), (18) can be rewritten as 

£ mSLX £ i min(m max ^) 

E E E E <-(-/ 2 ) 

£—0 m— — £ m' — — t m" — — min(m max ,^) 
i[m(a-7T /2)+m' j3+m" (7+tt /2)] 



TABLE I 

Algorithm complexity for one scale and one orientation. 



x e 



(21) 



where the symmetry relationship d £ mm , (— fl) = d^, m (/3) has 
been applied. (A similar factoring of the rotation and harmonic 
space representation has been independently performed in 
[48].) Steerable wavelets may have a low effective band limit 
m max <C 4ax, in which case the the inner summation in (21) 
may be truncated to m max . For general directional wavelets 
this is not the case and one must use m max = ^ ma x. 

Evaluating the harmonic formulation given by (21) directly 
would be no more efficient that approximating the initial inte- 
gral (8) using simple quadrature. However, (21) is represented 
in such a way that the presence of complex exponentials may 
be exploited such that FFTs may be applied to evaluate rapidly 
the three summations simultaneously. 

2) Fast implementation: Azimuthal rotations may be ap- 
plied with far less computational expense than polar rotations 
since they appear within complex exponentials in (21). Al- 
though the d-matrices can be evaluated reasonably quickly 
and reliably using recursion formulae, the basis for the fast 
implementation is to avoid these polar rotations as much as 
possible and use FFTs to evaluate rapidly all of the azimuthal 
rotations simultaneously. This is the motivation for factoring 
the rotation by (20) so that all Euler angles occur as azimuthal 
rotations. 

The discretisation of each Euler angle may in general be 
arbitrary. However, to exploit standard FFT routines uniform 
sampling is adopted, i.e. grid £2 is used (see section III- A). 
As mentioned, this discretisation is inefficient since it covers 
the SO (3) manifold exactly twice, nevertheless it enables the 
use of standard FFT routines which significantly increases the 
speed of the algorithm. Discretising (21) in this manner and 
interchanging the order of summation we obtain 

^-max ^max Trnax 

E E E 

-^max 

(*/2) r tm » 

^=max(|m| , \m'\ ,\m"\) 

v i[m(27tn a / N a -7r /2)+m'27rnB / Np+m" (27rn~ / N~+n /2)] 
X o 

Shifting the indices yields 

2^max 2£ max 2?n ma x 

e i2K(n a m/N a +nf3m' / Np+n^m" / N^) 



Algorithm 


Complexity 


Direct 


O(N^) 


Semi-fast 


G(N 3 log 2 N) 


Fast azimuthally symmetric 


G(N 3 ) 


Fast directional 


0(N 3 ) 



where 



rp _ i(m"-m)7r/2 

± m,m,m ff — c 



x E <L> m (*/2)dL> m »(*/2) 

^=max(|m|,|m'|,|m"|) 

is extended periodically. Note that the phase shift intro- 
duced by shifting the indices of the summation in (22), 
shifts the Tm^^i indices back. Making the associations 



1, N p = 2£ n 



1 and = 2m n 



1, one 



m=0 m'—O m"=0 



notices that (22) is the unnormalised 3 -dimensional inverse 
DFT of (23). Nyquist sampling in (a,/3, 7) is inherently 
ensured by the associations made with ^ max and m max . 

The CSWT may be performed rapidly in spherical harmonic 
space by constructing the T-matrix of (23) from spherical 
harmonic coefficients and precomputed d-matrices, followed 
by the application of an FFT to evaluate rapidly all three 
Euler angles of the discretised CSWT simultaneously. The 
complexity of the algorithm is dominated by the computation 
of the T-matrix. This involves performing a 1 -dimensional 
summation over a 3 -dimensional grid, hence the algorithm is 
of order 0(N 7 N 3 ). 

Additional benefits may be realised for real signals by ex- 
ploiting the resulting conjugate symmetry relationship. Mem- 
ory and computational requirements may be reduced by a 
further factor of two for real signals by noting the conjugate 
symmetry relationship T_ m _ m / _ m // = T m rn , m „. In our 
implementation we use the complex-to-real FFT routines of 
the FFTW 9 package, which are approximately twice as fast as 
the equivalent complex-to-complex routines. 

F. Comparison 

We summarise the computational complexities of the vari- 
ous CSWT algorithms for a single pair of scales and single 
orientation in Table I. The complexity of each algorithm 
scales with the number of dilations considered and, for those 
algorithms that facilitate a directional analysis {i.e. all but 
the fast azimuthally symmetric algorithm), with the number 
of orientations considered. The most general fast directional 
algorithm provides a saving of O (TV) over the direct case, 
where the number of pixels on the sphere and the harmonic 
band limit are related to N by G(N p - lyi ) = 0(£ max 2 ) = 
G(N 2 ). 

We implement all algorithms in Fortran 90, adopting the 
HEALPix tessellation of the sphere (which, incidentally, is 



(22) 9 http : / /www . f f tw . org/ 
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TABLE II 

Typical execution time (minutes :seconds) for each algorithm 
run on an Intel P4-M 3GHz laptop with 5 1 2MB of memory. 



Resolution 




Algorithm execution time 


^Vgide 


Npix 


Direct 


Fast isotropic 


Fast directional 


32 


12,288 


3:25.37 


0:00.06 


0:00.10 


64 


49,152 


54:31.75 


0:00.38 


0:00.74 


256 


786,432 




0:28.00 


0:52.55 


512 


3,257,292 




3:43.69 


7:57.75 


1024 


12,582,912 




28:23.85 


71:31.68 



the tessellation scheme of the WMAP CMB data [2]). Typical 
execution times for the algorithms are tabulated in Table II for 
a range of resolutions from 110 down to 3.4 arcminutes. The 
improvements provided by the fast algorithms are apparent. 
Indeed, it is not feasible to run the direct algorithm on data- 
sets with a resolution much greater than Api x c± 5 x 10 5 . For 
data-sets of practical size, such as the WMAP (iV P i x — 3 x 10 6 ) 
and forthcoming Planck (iV pix ~ 50 x 10 6 ) CMB data, the fast 
implementations of the CSWT are essential. 

The semi-fast algorithm is also implemented using the 
HEALPix tessellation. However, to perform the outer summa- 
tion (integration) continual interpolation followed by down- 
sampling is required on each iso-latitude ring to essentially 
resample the data on an ECP tessellation. This increases 
the execution time of the implementation of the semi-fast 
algorithm on the HEALPix grid to an extent that the semi-fast 
algorithm provides little advantage over the direct algorithm. 
To appreciate the advantages of the semi-fast approach it 
must be implemented on an ECP tessellation, hence we do 
not tabulate the execution times for our implementation of 
this algorithm on the HEALPix grid as it provides an unfair 
comparison. 

It is also important to note that although complexity scales 
with the number of dilations and orientations considered, 
execution time does not for the fast algorithms. Execution 
time does scale in this manner for the direct algorithm. 
There are a number of additional overheads associated with 
the fast algorithms, such as computing spherical harmonic 
coefficients and <i-matrices. Consequently, the fast algorithms 
provide additional execution time savings that are not directly 
apparent in Table II. For example, the execution time of the 
fast azimuthally symmetric and directional algorithms for 10 
dilations at a resolution of Api x ^8x 10 5 (A^de = 256) are 
3:08.20 and 7:06.83 (minutes: seconds) respectively, which is 
considerably faster than ten times the execution time of one 
dilation. 

IV. Application 

We demonstrate in this section the application of our CSWT 
implementation to binary Earth data. In Fig. 3 the Earth 
data and the corresponding spherical wavelet coefficients are 
shown. We use the SBW defined in section II-D to perform 
the analysis. This is a steerable wavelet [27], however our 
implementation is in general valid for any directional wavelet. 
Notice how the wavelet coefficient maps corresponding to 
different oriented wavelets pick out corresponding oriented 



structure in the data. As the dilation scale is increased, the 
scale of the features extracted also increases accordingly. 

The ability to probe oriented structure in data defined 
on a sphere is of important practical use. Certain physical 
processes may be localised on the sphere in scale, position 
and orientation (e.g. signatures of cosmic strings in the CMB 
[49] or correlations induced in the CMB by the nearby galaxy 
distribution [50]). Thus, analysing the statistical properties 
of spherical wavelet coefficients individually for a range of 
scales, positions and orientations may allow one to detect 
such effects with greater significance. Indeed, using a direc- 
tional spherical wavelet analysis we have made very strong 
detections of non-Gaussianity in the CMB [36], [37] and the 
strongest detection made to date of the ISW effect [38]. 

V. Concluding remarks 

The extension of Euclidean wavelet analysis to the sphere 
has been described in the framework presented by [27], where 
the stereographic projection is used to relate the spherical 
and Euclidean constructions. We extend the concept of the 
spherical dilation presented by [27] to anisotropic dilations. 
Although anisotropic dilations are of practical use, the result- 
ing basis one projects onto does not classify as a wavelet basis 
since perfect reconstruction is not possible. 

Current and forthcoming data-sets on the sphere, of the 
CMB for example, are of considerable size as higher reso- 
lutions are necessary to test new physics. Consequently, we 
present fast algorithms to implement the CSWT as an analysis 
without such algorithms is not computationally feasible. A 
range of algorithms are described, from the direct quadra- 
ture approximation, to the semi-fast equi-angular algorithm 
where one rotation is performed in Fourier space, to the 
fast azimuthally symmetric and directional algorithms posed 
purely in spherical harmonic space. Posing the CSWT purely 
in harmonic space naturally ensures the resulting algorithms 
are tessellation independent. The most general fast directional 
algorithm provides a saving of O(^fN^) = 0(£ max ) over 
the direct implementation and may be performed down to a 
few arcminutes even with limited computational resources. 

Data is observed on a sphere in a range of applications. In 
many of these cases the ability to perform a wavelet analysis 
would be useful. For example, spherical wavelets may be used 
to probe the CMB for deviations from the standard assumption 
of Gaussianity, or to search for compact objects embedded 
in the CMB, such as cosmic strings, a predicted but as yet 
unobserved phenomenon. The extension of wavelet analysis 
to the sphere enables one to probe new physics in many 
areas, and is facilitated on large practical data-sets by our fast 
directional CSWT algorithm. 
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(a) Binary Earth map 




(b) Spherical butterfly wavelet coefficients for 
scale a = b = 0.03 and orientation 7 = 0° 




(e) Spherical butterfly wavelet coefficients for 
scale a = b = 0.12 and orientation 7 = 144° 

Fig. 3. Binary Earth data and corresponding spherical butterfly wavelet 
coefficients. The wavelet coefficients for combinations of only two scales 
and two orientations are shown. Notice how the wavelet coefficient maps 
corresponding to different oriented wavelets pick out corresponding oriented 
strnrtiire in the. data As the. dilation srale. is increased the. srale. of the. features 



use of the YAWTb Matlab toolbox for the binary Earth data 
defined on the sphere. 
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